R Notebook

library(TSA)
library(tseries)
library(astsa)
library(imputeTS)
library(tsoutliers)
library(xts)

Interpolating Data

Here we interpolate our missing data with a linear model.

terror2 <- read.csv("input/og_num_casualities_greater_than_10.csv")
terror3 <- na.interpolation(terror2$num.attacks.with.kill.thresh, option="linear")
plot(as.xts(ts(terror3, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (w/ Linear Imputed Data)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", ylim=c(0, 225), col="red")
pdf("image/og_ts.pdf")

lines(as.xts(ts(terror2$num.attacks.with.kill.thresh, frequency = 12, start=1970)), col="black", lwd=1.2)
dev.off()

quartz_off_screen 
                2 

Removing outliers

outlier_terror3 <- tso(ts(terror3), types = c("TC", "AO", "IO"))

stopped when 'maxit.oloop = 4' was reached

plot(outlier_terror3)
#plot outlier effects
pdf("image/outlier_effects.pdf")

plot(as.xts(ts(outlier_terror3$effects, frequency = 12, start=1970)), main = "Outlier Effects", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="red")
dev.off()

quartz_off_screen 
                2 

#Plot outlier time series
xts.terror3 <- as.xts(ts(terror3, frequency = 12, start=1970))
plot(as.xts(ts(terror3, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Outliers Removed)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="lightgray", ylim=c(0, 225))

lines(as.xts(ts(outlier_terror3$yadj, frequency = 12, start=1970)), col="blue")

points(xts.terror3[427], col="red",pch=19, cex=1)

points(xts.terror3[516], col="red",pch=19, cex=1)

points(xts.terror3[521], col="red",pch=19, cex=1)

points(xts.terror3[523], col="red",pch=19, cex=1)

points(xts.terror3[547], col="red",pch=19, cex=1)
pdf("image/outlier_comparison.pdf")

points(xts.terror3[556], col="red",pch=19, cex=1)
dev.off()

quartz_off_screen 
                2 

Making, Training, and Validation set

terror4 <- outlier_terror3$yadj
#terror3 <- na.kalman(terror2$num.attacks, model="auto.arima")
cuttoff.index <- length(terror4) - 24 #floor(0.1 * length(terror3))
cuttoff.index2 <- length(terror4) - 12
terror4.valid <- terror4[(cuttoff.index+1) :cuttoff.index2]
terror4.testing <- terror4[(cuttoff.index2 + 1): length(terror4)]
terror4 <- terror4[1: cuttoff.index]
#plot(as.xts(ts(terror4, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Training Set)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
#plot(as.xts(ts(terror4.valid, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Validation Set)", major.format = "%Y-%m", grid.col="white", lwd=1)

Chasing Stationarity

#log_terror4 <- log(outlier_terror3$yadj)
adf.test(terror4, k=1)

    Augmented Dickey-Fuller Test

data:  terror4
Dickey-Fuller = -2.4189, Lag order = 1, p-value = 0.401
alternative hypothesis: stationary

adf.test(diff(terror4), k=1)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  diff(terror4)
Dickey-Fuller = -24.444, Lag order = 1, p-value = 0.01
alternative hypothesis: stationary

adf.test(diff(diff(terror4)), k=1)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  diff(diff(terror4))
Dickey-Fuller = -31.986, Lag order = 1, p-value = 0.01
alternative hypothesis: stationary

pdf("image/first_diff.pdf")
plot(as.xts(ts(diff(terror4), frequency = 12, start=1970)), main = "Number of Terrorist Attacks (First Diff)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
dev.off()

null device 
          1 

pdf("image/second_diff.pdf")

plot(as.xts(ts(diff(diff(terror4)), frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Second Diff)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
dev.off()

null device 
          1 

#ts.plot(diff(terror4))
#ts.plot(diff(diff(terror4)))

pdf("image/acf_og.pdf")
acf(terror4, main="ACF of Training Data")
dev.off()

null device 
          1 

pdf("image/pacf_og.pdf")

pacf(terror4, main="PACF of Training Data")
dev.off()

null device 
          1 

pdf("image/acf_first_diff.pdf")

acf(diff(terror4), main="ACF of First Diff Training Data")
dev.off()

null device 
          1 

pdf("image/pacf_first_diff.pdf")

pacf(diff(terror4), main="PACF of First Diff Training Data")
dev.off()

null device 
          1 

pdf("image/acf_second_diff.pdf")

acf(diff(terror4), main="ACF of Second Diff Training Data")
dev.off()

null device 
          1 

pdf("image/pacf_second_diff.pdf")

pacf(diff(terror4), main="PACF of Second Diff Training Data")
dev.off()

null device 
          1 

Periodogram; Figuring out Seasonality

m = floor(sqrt(length(diff(terror4))))
#pdf("image/raw_periodogram.pdf")
spec.pgram(diff(terror4), log="no", main="Raw Periodogram", cex.main=1.5)

#dev.off()
#pdf("image/smooth_tapered_periodogram.pdf")
spec.pgram(diff(terror4), kernel('daniell', m), log="no", taper=0.1, main="Smoothed and Tapered Periodogram")

#dev.off()
#mvspec(diff(log_terror4),  kernel('daniell', m), log="no")

Finding which model to use

eacf(diff(terror4))

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o x x o o o o o o o  o  o  x 
1 x o o x o o o o o o o  o  o  o 
2 x x x x o o o o o o o  o  x  o 
3 o x x x o o o o o o o  o  x  o 
4 x o o o x o x o o o o  o  o  o 
5 x x o o o o o o o o o  o  o  o 
6 x x x o o o o o o o o  o  o  o 
7 x x x x o o o o o o o  o  o  o 

eacf(diff(diff(terror4)))

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x x x x o o o o o o o  o  o  x 
1 x x x x o o o o o o o  o  o  x 
2 x x x x o x o o o o o  o  o  o 
3 x x o o o o o o o o o  o  o  o 
4 x x o o o o x x o o o  o  o  o 
5 x x o o o o o o o o o  o  o  o 
6 x o o x o o o o o o o  o  o  o 
7 x x o x x o o o o o o  o  o  o 

#sarima(terror4, 0, 1, 1)
sarima(terror4, 0, 1, 1, 1, 0, 1, 4)

initial  value 2.346947 
iter   2 value 2.220389
iter   3 value 2.211951
iter   4 value 2.193176
iter   5 value 2.187658
iter   6 value 2.185821
iter   7 value 2.185525
iter   8 value 2.185512
iter   9 value 2.185508
iter  10 value 2.185484
iter  11 value 2.185474
iter  12 value 2.185467
iter  13 value 2.185459
iter  14 value 2.185419
iter  15 value 2.185268
iter  16 value 2.184298
iter  17 value 2.183991
iter  18 value 2.182728
iter  19 value 2.181869
iter  20 value 2.181763
iter  21 value 2.181685
iter  22 value 2.181580
iter  23 value 2.181485
iter  24 value 2.181480
iter  25 value 2.181480
iter  25 value 2.181480
iter  25 value 2.181480
final  value 2.181480 
converged
initial  value 2.178465 
iter   2 value 2.178461
iter   3 value 2.178459
iter   4 value 2.178459
iter   5 value 2.178458
iter   6 value 2.178457
iter   7 value 2.178457
iter   8 value 2.178457
iter   8 value 2.178457
iter   8 value 2.178457
final  value 2.178457 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
    reltol = tol))

Coefficients:
          ma1     sar1    sma1  constant
      -0.6128  -0.6221  0.7609    0.2964
s.e.   0.0351   0.1214  0.1003    0.1602

sigma^2 estimated as 77.91:  log likelihood = -1939,  aic = 3887.99

$degrees_of_freedom
[1] 535

$ttable
         Estimate     SE  t.value p.value
ma1       -0.6128 0.0351 -17.4617  0.0000
sar1      -0.6221 0.1214  -5.1236  0.0000
sma1       0.7609 0.1003   7.5870  0.0000
constant   0.2964 0.1602   1.8496  0.0649

$AIC
[1] 5.370387

$AICc
[1] 5.374299

$BIC
[1] 4.402176

sarima(terror4, 0, 1, 1, 1, 1, 1, 4)

initial  value 2.599103 
iter   2 value 2.332712
iter   3 value 2.270429
iter   4 value 2.237537
iter   5 value 2.222895
iter   6 value 2.215065
iter   7 value 2.205793
iter   8 value 2.204129
iter   9 value 2.203212
iter  10 value 2.203071
iter  11 value 2.202956
iter  12 value 2.202950
iter  13 value 2.202949
iter  13 value 2.202949
iter  13 value 2.202949
final  value 2.202949 
converged
initial  value 2.204741 
iter   2 value 2.203758
iter   3 value 2.203347
iter   4 value 2.203343
iter   5 value 2.203343
iter   6 value 2.203342
iter   6 value 2.203342
iter   6 value 2.203342
final  value 2.203342 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc, 
    REPORT = 1, reltol = tol))

Coefficients:
          ma1    sar1     sma1
      -0.6121  0.1591  -1.0000
s.e.   0.0375  0.0434   0.0362

sigma^2 estimated as 79.13:  log likelihood = -1937.92,  aic = 3883.84

$degrees_of_freedom
[1] 532

$ttable
     Estimate     SE  t.value p.value
ma1   -0.6121 0.0375 -16.3203   0e+00
sar1   0.1591 0.0434   3.6675   3e-04
sma1  -1.0000 0.0362 -27.6028   0e+00

$AIC
[1] 5.382182

$AICc
[1] 5.386024

$BIC
[1] 4.406024

sarima(terror4, 0, 1, 1, 1, 1, 2, 4)

initial  value 2.599103 
iter   2 value 2.300969
iter   3 value 2.258534
iter   4 value 2.236661
iter   5 value 2.205300
iter   6 value 2.201253
iter   7 value 2.198386
iter   8 value 2.197243
iter   9 value 2.196562
iter  10 value 2.196267
iter  11 value 2.196213
iter  12 value 2.196133
iter  13 value 2.195850
iter  14 value 2.194818
iter  15 value 2.194074
iter  16 value 2.193350
iter  17 value 2.192729
iter  18 value 2.192231
iter  19 value 2.192117
iter  20 value 2.192100
iter  21 value 2.192094
iter  22 value 2.192032
iter  23 value 2.191993
iter  24 value 2.191956
iter  25 value 2.191945
iter  26 value 2.191944
iter  26 value 2.191944
iter  26 value 2.191944
final  value 2.191944 
converged
initial  value 2.199972 
iter   2 value 2.199945
iter   3 value 2.199936
iter   4 value 2.199929
iter   5 value 2.199892
iter   6 value 2.199876
iter   7 value 2.199872
iter   8 value 2.199872
iter   8 value 2.199872
iter   8 value 2.199872
final  value 2.199872 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc, 
    REPORT = 1, reltol = tol))

Coefficients:
          ma1     sar1     sma1     sma2
      -0.6168  -0.6197  -0.2244  -0.7477
s.e.   0.0381   0.1241   0.1090   0.0999

sigma^2 estimated as 79.3:  log likelihood = -1936.06,  aic = 3882.13

$degrees_of_freedom
[1] 531

$ttable
     Estimate     SE  t.value p.value
ma1   -0.6168 0.0381 -16.2046    0.00
sar1  -0.6197 0.1241  -4.9951    0.00
sma1  -0.2244 0.1090  -2.0593    0.04
sma2  -0.7477 0.0999  -7.4861    0.00

$AIC
[1] 5.388087

$AICc
[1] 5.391999

$BIC
[1] 4.419876

sarima(terror4, 1, 1, 2)

initial  value 2.344191 
iter   2 value 2.241534
iter   3 value 2.195978
iter   4 value 2.194817
iter   5 value 2.194230
iter   6 value 2.194208
iter   7 value 2.194208
iter   8 value 2.194207
iter   9 value 2.194207
iter  10 value 2.194205
iter  11 value 2.194205
iter  12 value 2.194204
iter  13 value 2.194204
iter  14 value 2.194204
iter  15 value 2.194201
iter  16 value 2.194195
iter  17 value 2.194183
iter  18 value 2.194180
iter  19 value 2.194178
iter  20 value 2.194178
iter  21 value 2.194178
iter  21 value 2.194178
iter  21 value 2.194178
final  value 2.194178 
converged
initial  value 2.193663 
iter   2 value 2.193663
iter   3 value 2.193663
iter   3 value 2.193663
iter   3 value 2.193663
final  value 2.193663 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
    reltol = tol))

Coefficients:
          ar1      ma1      ma2  constant
      -0.0950  -0.4967  -0.0688    0.2960
s.e.   0.7203   0.7171   0.4282    0.1536

sigma^2 estimated as 80.36:  log likelihood = -1947.19,  aic = 3904.38

$degrees_of_freedom
[1] 535

$ttable
         Estimate     SE t.value p.value
ar1       -0.0950 0.7203 -0.1320  0.8951
ma1       -0.4967 0.7171 -0.6926  0.4888
ma2       -0.0688 0.4282 -0.1608  0.8723
constant   0.2960 0.1536  1.9265  0.0546

$AIC
[1] 5.401315

$AICc
[1] 5.405227

$BIC
[1] 4.433105

sarima(terror4, 1, 1, 2, 1, 0, 1, 4)

initial  value 2.347868 
iter   2 value 2.258620
iter   3 value 2.188794
iter   4 value 2.186745
iter   5 value 2.185962
iter   6 value 2.185836
iter   7 value 2.185771
iter   8 value 2.185767
iter   9 value 2.185742
iter  10 value 2.185679
iter  11 value 2.185479
iter  12 value 2.185171
iter  13 value 2.184101
iter  14 value 2.183062
iter  15 value 2.182366
iter  16 value 2.182283
iter  17 value 2.182258
iter  18 value 2.182224
iter  19 value 2.182152
iter  20 value 2.182073
iter  21 value 2.182032
iter  22 value 2.181984
iter  23 value 2.181979
iter  24 value 2.181979
iter  24 value 2.181979
iter  24 value 2.181979
final  value 2.181979 
converged
initial  value 2.178050 
iter   2 value 2.178045
iter   3 value 2.178044
iter   4 value 2.178044
iter   5 value 2.178044
iter   6 value 2.178043
iter   7 value 2.178043
iter   8 value 2.178043
iter   8 value 2.178043
iter   8 value 2.178043
final  value 2.178043 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
    reltol = tol))

Coefficients:
         ar1      ma1     ma2     sar1    sma1  constant
      0.1981  -0.7935  0.0937  -0.6070  0.7515    0.2966
s.e.  0.8900   0.8938  0.5535   0.1345  0.1104    0.1555

sigma^2 estimated as 77.85:  log likelihood = -1938.77,  aic = 3891.55

$degrees_of_freedom
[1] 533

$ttable
         Estimate     SE t.value p.value
ar1        0.1981 0.8900  0.2226  0.8239
ma1       -0.7935 0.8938 -0.8878  0.3750
ma2        0.0937 0.5535  0.1692  0.8657
sar1      -0.6070 0.1345 -4.5117  0.0000
sma1       0.7515 0.1104  6.8085  0.0000
constant   0.2966 0.1555  1.9074  0.0570

$AIC
[1] 5.376974

$AICc
[1] 5.381068

$BIC
[1] 4.424658

sarima(terror4, 1, 1, 2, 1, 1, 1, 4)

initial  value 2.600040 
iter   2 value 2.356718
iter   3 value 2.277980
iter   4 value 2.233407
iter   5 value 2.223552
iter   6 value 2.212201
iter   7 value 2.205827
iter   8 value 2.203297
iter   9 value 2.201565
iter  10 value 2.200773
iter  11 value 2.200719
iter  12 value 2.200705
iter  13 value 2.200661
iter  14 value 2.200544
iter  15 value 2.200182
iter  16 value 2.200058
iter  17 value 2.199536
iter  18 value 2.198925
iter  19 value 2.198604
iter  20 value 2.198391
iter  21 value 2.198352
iter  22 value 2.198279
iter  23 value 2.198233
iter  24 value 2.197691
iter  25 value 2.195332
iter  26 value 2.195248
iter  27 value 2.193232
iter  28 value 2.191968
iter  29 value 2.190157
iter  30 value 2.190012
iter  31 value 2.189751
iter  32 value 2.189563
iter  33 value 2.189216
iter  34 value 2.188867
iter  35 value 2.188828
iter  36 value 2.188810
iter  37 value 2.188802
iter  37 value 2.188802
iter  37 value 2.188802
final  value 2.188802 
converged
initial  value 2.202122 
iter   2 value 2.202104
iter   3 value 2.202047
iter   4 value 2.202036
iter   5 value 2.202000
iter   6 value 2.201960
iter   7 value 2.201834
iter   8 value 2.201715
iter   9 value 2.201577
iter  10 value 2.201560
iter  11 value 2.201548
iter  12 value 2.201543
iter  12 value 2.201543
final  value 2.201543 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc, 
    REPORT = 1, reltol = tol))

Coefficients:
         ar1      ma1     ma2    sar1     sma1
      0.4461  -1.0292  0.2206  0.1873  -0.9933
s.e.  0.3240   0.3308  0.2148  0.0489   0.0407

sigma^2 estimated as 79.29:  log likelihood = -1936.96,  aic = 3885.92

$degrees_of_freedom
[1] 530

$ttable
     Estimate     SE  t.value p.value
ar1    0.4461 0.3240   1.3767  0.1692
ma1   -1.0292 0.3308  -3.1116  0.0020
ma2    0.2206 0.2148   1.0270  0.3049
sar1   0.1873 0.0489   3.8268  0.0001
sma1  -0.9933 0.0407 -24.3931  0.0000

$AIC
[1] 5.39168

$AICc
[1] 5.395676

$BIC
[1] 4.431417

sarima(terror4, 1, 1, 2, 1, 1, 2, 4)

initial  value 2.600040 
iter   2 value 2.325435
iter   3 value 2.237416
iter   4 value 2.214251
iter   5 value 2.202589
iter   6 value 2.198408
iter   7 value 2.198339
iter   8 value 2.197318
iter   9 value 2.196580
iter  10 value 2.196494
iter  11 value 2.196379
iter  12 value 2.196322
iter  13 value 2.195990
iter  14 value 2.195374
iter  15 value 2.195370
iter  16 value 2.195119
iter  17 value 2.194986
iter  18 value 2.194830
iter  19 value 2.194773
iter  20 value 2.194505
iter  21 value 2.194160
iter  22 value 2.193795
iter  23 value 2.193655
iter  24 value 2.193455
iter  25 value 2.193431
iter  26 value 2.193425
iter  27 value 2.193424
iter  27 value 2.193424
final  value 2.193424 
converged
initial  value 2.199628 
iter   2 value 2.199619
iter   3 value 2.199609
iter   4 value 2.199604
iter   5 value 2.199600
iter   6 value 2.199580
iter   7 value 2.199565
iter   8 value 2.199563
iter   9 value 2.199562
iter  10 value 2.199561
iter  11 value 2.199557
iter  12 value 2.199550
iter  13 value 2.199534
iter  14 value 2.199505
iter  15 value 2.199481
iter  16 value 2.199471
iter  17 value 2.199467
iter  18 value 2.199467
iter  19 value 2.199467
iter  20 value 2.199464
iter  21 value 2.199458
iter  22 value 2.199449
iter  23 value 2.199441
iter  24 value 2.199438
iter  25 value 2.199436
iter  26 value 2.199433
iter  27 value 2.199426
iter  28 value 2.199422
iter  29 value 2.199420
iter  30 value 2.199419
iter  31 value 2.199417
iter  32 value 2.199414
iter  33 value 2.199410
iter  34 value 2.199407
iter  35 value 2.199407
iter  36 value 2.199407
iter  37 value 2.199406
iter  38 value 2.199406
iter  39 value 2.199404
iter  40 value 2.199400
iter  40 value 2.199400
final  value 2.199400 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc, 
    REPORT = 1, reltol = tol))

Coefficients:
         ar1      ma1     ma2     sar1     sma1     sma2
      0.0787  -0.6766  0.0177  -0.6062  -0.2298  -0.7358
s.e.  0.9149   0.9129  0.5568   0.1362   0.1180   0.1096

sigma^2 estimated as 79.34:  log likelihood = -1935.81,  aic = 3885.62

$degrees_of_freedom
[1] 529

$ttable
     Estimate     SE t.value p.value
ar1    0.0787 0.9149  0.0860  0.9315
ma1   -0.6766 0.9129 -0.7412  0.4589
ma2    0.0177 0.5568  0.0318  0.9746
sar1  -0.6062 0.1362 -4.4501  0.0000
sma1  -0.2298 0.1180 -1.9474  0.0520
sma2  -0.7358 0.1096 -6.7148  0.0000

$AIC
[1] 5.395983

$AICc
[1] 5.400077

$BIC
[1] 4.443667

sarima(terror4, 1, 1, 1)

initial  value 2.344191 
iter   2 value 2.241424
iter   3 value 2.210360
iter   4 value 2.203987
iter   5 value 2.203217
iter   6 value 2.196675
iter   7 value 2.194860
iter   8 value 2.194317
iter   9 value 2.194264
iter  10 value 2.194223
iter  11 value 2.194209
iter  12 value 2.194204
iter  13 value 2.194204
iter  13 value 2.194204
iter  13 value 2.194204
final  value 2.194204 
converged
initial  value 2.193680 
iter   2 value 2.193680
iter   3 value 2.193679
iter   3 value 2.193679
iter   3 value 2.193679
final  value 2.193679 
converged
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
    Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
    reltol = tol))

Coefficients:
         ar1      ma1  constant
      0.0175  -0.6101    0.2960
s.e.  0.0678   0.0517    0.1537

sigma^2 estimated as 80.36:  log likelihood = -1947.2,  aic = 3902.4

$degrees_of_freedom
[1] 536

$ttable
         Estimate     SE  t.value p.value
ar1        0.0175 0.0678   0.2580  0.7965
ma1       -0.6101 0.0517 -11.7941  0.0000
constant   0.2960 0.1537   1.9258  0.0547

$AIC
[1] 5.397646

$AICc
[1] 5.401488

$BIC
[1] 4.421488

pdf("image/best_model.pdf")
sarima(terror4, 0, 1, 1, 1, 1, 1, 4)
dev.off()

MSE calculations

#eacf(diff(diff(log_terror4)))
#terror5 <- terror4
#total_error <- 0
#for (i in 1: (length(terror4.valid) - 11))
#{
#  actual <- terror4.valid[i : i + 11]
#  predicted <- sarima.for(terror5, 12, 0, 1, 1, 0, 0, 0, 0) $pred
#  total_error <- total_error + sum((actual - predicted)^2)
#  terror5 <- c(terror5, terror4.valid[i]) 
#}
predicted <- sarima.for(terror4, 12, 0, 1, 1, 0, 0, 0, 0)$pred

mse <- sum((predicted - terror4.valid)^2)
mse

[1] 4248.281

Predicting the future using our best model

val <- sarima.for(c(terror4, terror4.valid), 12, 0, 1, 1, 1, 1, 1, 4)
pred <-val$pred
err  <-val$se
total <- c(terror4, terror4.valid, terror4.testing)
par(cex.main = 2)

plot(as.xts(ts(total, frequency = 12, start=1970))[492:length(total)], main = "Number of Terrorist Attacks (Prediction)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="lightgray", pch="1", ylim=c(0, 225))

points(as.xts(ts(total, frequency = 12, start=1970)),col="lightgray",pch="o")

lines(as.xts(ts(c(terror4, terror4.valid), frequency = 12, start=1970)),col="black")

points(as.xts(ts(c(terror4, terror4.valid), frequency = 12, start=1970)),col="black",pch="o")

lines(as.xts(ts(pred, frequency = 12, start=2016)),col="blue")

lines(as.xts(ts(pred + err, frequency = 12, start=2016)),col="blue", lty="dashed")

lines(as.xts(ts(pred - err, frequency = 12, start=2016)),col="blue", lty="dashed")

lines(as.xts(ts(pred + 2*err, frequency = 12, start=2016)),col="blue", lty="dotted")
pdf("image/prediction_on_testing.pdf")

lines(as.xts(ts(pred - 2*err, frequency = 12, start=2016)),col="blue", lty="dotted")
dev.off()

quartz_off_screen 
                2 

mse <- sum((pred - terror4.testing)^2)
mse

[1] 1784.668